from numpy import *
from numpy.linalg import *
from pylab import *
import matplotlib.pyplot as mp

#fonctions exemple
def fon(t, y):
	return y/(1+t**2)

def fon2(t, y):
	return array([-y[1],y[0]])


                               ############################
                               # FONCTIONS DE LA PARTIE 1 #
                               ###########################

#methodes pour un pas
def step_euler(f,t0,y0,h):
	t1=t0+h
	y1=y0+h*f(t0,y0)
	return (t1, y1)

def step_pt_milieu(f,t0,y0,h):
	t1=t0+h
	y1=y0+h*f(t0+h/2., y0+h/2.*f(t0,y0))
	return (t1, y1)

def step_heun(f,t0,y0,h):
	t1=t0+h
	y1=y0+h/2.*(f(t0,y0)+f(t0+h,y0+h*f(t0,y0)))
	return (t1, y1)

def step_rk4(f,t0,y0,h):
        k1=f(t0,y0)
        k2=f(t0+h/2, y0+h/2*k1)
        k3=f(t0+h/2, y0+h/2*k2)
        k4=f(t0+h, y0+h*k3)
	t1=t0+h
        y1=y0+h/6*(k1+2.*k2+2.*k3+k4)
	return (t1, y1)


#fonction generale pour un nombre de pas fixe
def meth_n_step(f,  #fonction du type f(t, y)
		t0, #debut de l'intervalle considere (int/float)
		y0, #cond init (int/float ou array/liste si ED de dim>1)
		h,  #pas (int/float)
		N,  #nombre de pas (int/float)
		meth="rk4", #methode ("euler", "pt_milieu", "heun" ou "rk4")
		aff="non"   #afficher courbe ou retourner resultats ("oui"/rien)
		):

	#changement des types des variables
	typ_y0=type(y0)
	t0=float(t0)
	h=float(h)
	bo=False
	if typ_y0==ndarray or typ_y0==list:
		bo=True
		if typ_y0==list:
			y0=array(y0)
		for i in range(0, len(y0)):
			y0[i]=float(y0[i])

	X=[t0]
	Y=[y0]

	#choix de la methode
	if meth=="euler":
		etap=step_euler
	elif meth=="pt_milieu":
		etap=step_pt_milieu
	elif meth=="heun":
		etap=step_heun
	elif meth=="rk4":
		etap=step_rk4

	#boucle des iterations de la methode
	while N!=0:
		A=etap(f, t0, y0, h)

		#changement des types si necessaire
		if bo:
			if type(A[0])==list:
				A[0]=array(A[0])
			if type(A[1])==list:
				A[1]=array(A[1])

		X.append(A[0])
		Y.append(A[1])
		t0=A[0]
		y0=A[1]
		N=N-1

	#affichage ou retour des resultats
	if aff=="oui":
		mp.plot(X,Y)
		mp.savefig("figure")
		mp.show()
	else:
		return (X,Y)


#fonction generale pour une precision demandee
def meth_epsilon(f,  #fonction du type f(t, y)
		 y0, #condition initiale (type array ou liste si ED de dim>1)
		 t0, #debut de l'intervalle considere
		 tf, #fin de l'intervalle considere
		 eps=0.01,   #precision
		 meth="rk4", #methode ("euler", "pt_milieu", "heun" ou "rk4")
		 aff="oui"   #afficher courbe ou retourner resultats
		 ):
	bo=True
	typ=type(y0)
	bo2=typ==ndarray or typ==list
	if bo2:
		eps=array(len(y0)*[eps])
	else:
		y0=float(y0)
	t0=float(t0)
	tf=float(tf)
	N=10
	t=tf-t0
	while bo:
		h=t/N
		A=meth_n_step(f, t0, y0, h, N, meth, "non")
		A1=array(A[1])
		B=meth_n_step(f, t0, y0, h/2., N*2, meth, "non")
		B1=array(B[1])
		C=[]
		for i in range(0, len(A1)):
			C.append(B1[2*i])
		C=array(C)
		D=abs(C-A1)
		if bo2:
			m=array(len(y0)*[0.])
		else:
			m=0.
		for i in range(1, len(D)):
			m=m+D[i]
		m=m/N
		if (m<eps).all():
			bo=False
		N=N*2
	if aff=="oui":
		print "Nb de pas =",N
		print "precision =",m
		courbe = mp.plot(B[0],B[1])
		mp.legend(courbe ,("X","Y","Z"))
		xlabel("t")
		ylabel("y(t)")
		mp.savefig("figure")
		mp.show()
	else:
		return B

	###########################################partie3####################################################################
###########################################################################################################################################

B=linspace(0.1,0.2,20)
I=array([1,0,0])

def definir_sys(b):
    a=10
    c=200
    d=50
    
    def Zhabotinsky (t,N):
        return array([a+b*N[2]-N[0]-N[0]*N[1]*N[1],
                      c*(N[0]+N[0]*N[1]*N[1]-N[1]),
                      d*(N[1]-N[2])])
    return Zhabotinsky
    
def resoudre(I,b,eps,debut,fin):
        return meth_epsilon(definir_sys(b),I,debut,fin,eps,"rk4","non")



def partie_periodique(f,begin,x):
    e=10**-5
    #e est une tolerence pour eviter de detecter trop de pics
    t=begin
    fin=size(x)
        
    while t<fin-1 and f[t][0] - f[t+1][0]<e:
        t=t+1
    return (t)

def periode(I,b,eps,debut,end):
    (x,f)=resoudre(I,b,eps,debut,end)
    fin=size(x)
    pic=[]
    periodes=0
    securite=0
    trouve=1
    d=partie_periodique(f,0,x)
    #debut recoit la premiere abscisse de la partie periodique
    while d < fin-2  and securite < fin:
        
        pic=pic+[f[d][0]]
        d=partie_periodique(f,d+1,x)
        periodes=periodes+1
        securite=securite+1
    return (periodes,pic)

def bifurcation (I,B,eps,debut,fin) :
    for i in range(np.size(B)):
        X=periode(I,B[i],eps,debut,fin)[1]
        mp.plot([B[i]],[X],'.')
        print(i)
    xlabel("Coefficient b")
    ylabel("pics en ordonnee")
    mp.savefig("bifurcation")
    mp.show()
        
        




